set more off

clear matrix
clear mata
clear
set mem 10000m
set maxvar 32767
set matsize 11000

capture log close



global home "C:\Florida projects\dissertation\schools wind pollution"
global output "C:\Users\Dropbox\Research on Florida Wind Patters\School wind pollution\output\whole sample results"
global samples "$home\samples"


log using "$output/Figures112618_$S_DATE.log", replace

use "$samples\redceclptestwk_32219.dta", clear // BASED ON setuprdce_origtestweek_32219.do

if `David' == 1{

gen placebo_any_mjrhwy = runiform() if down60_any_mjrhwy_4~=.
gen absent_days = runiform() if down60_any_mjrhwy_4~=.
gen rabsent  = runiform() if down60_any_mjrhwy_4~=.
gen down60_any_mjrhwy_4_test = runiform() if down60_any_mjrhwy_4~=.
gen dwintensity4h_test = runiform() if down60_any_mjrhwy_4~=.
gen dwind_any_win4_test = runiform() if down60_any_mjrhwy_4~=.
}



if `Jenni'==1 {
	gen schnam="check"
	gen down60_any_mjrhwy_4=downwindmorehalf5
	gen dwind_any_win4_mjrhwy=downwind1
	gen percenthispanic=runiform()
	replace FRL=FRL/100
	bysort id (year): gen mover = school -school[_n-1]
		replace mover=1 if mover!=0
		bysort id (year): replace mover=0 if _n==1
	gen distmj1=frl
	gen rdcnt_win10_mjrhwy =frl
	gen behaveincident = (nincident<200)
	gen  graderep = behaveincident	
	gen AADT1_mjrhwy=runiform()*100000
	gen closetohwydwind60_4=black
	gen mitohwy=runiform()
	gen dwintensity4h=runiform()
	}


if `Jennihome'==1 {
	gen down60_any_mjrhwy_4=downwindmorehalf5
	gen dwind_any_win4_mjrhwy=runiform()
	gen percenthispanic=runiform()
	replace FRL=FRL/100
	bysort id (year): gen mover = school -school[_n-1]
		replace mover=1 if mover!=0
		bysort id (year): replace mover=0 if _n==1
	gen distmj1=frl
	gen rdcnt_win10_mjrhwy =frl
	gen zip=DistrictNumber
		destring zip, replace
	gen behaveincident = (nincident<200)
	gen  graderep = behaveincident
	gen rabsent=runiform()
	gen changein6or9id=1
	}

***********models*********************

global model3 momedbyschool momblackbyschool mommarriedbyschool percenthispanic frl  TeachDegree size stability mover distmj* rdcnt_win10_mjrhwy


***************************************


drop  if strpos(schnam, "VIRTUAL")  | strpos(schnam, "ONLINE") | strpos(schnam, "ADULT") | strpos(schnam, "AMIKIDS") ///
 | strpos(schnam, "AMI KIDS") | strpos(schnam, "JAIL") | strpos(schnam, "JUVENILE") | strpos(schnam, "TRANSITION") ///
 | strpos(schnam, "HALFWAY") | strpos(schnam, "LIFE SKILLS") | strpos(schnam, "DROPOUT") | strpos(schnam, "HOSPITAL") ///
  | strpos(schnam, "PARENT")  | strpos(schnam, "DETENTION") | strpos(schnam, "SUMMER") | strpos(schnam, "ADDICTIONS") ///
  | strpos(schnam, "PACE") | strpos(schnam, "SEAL CENTER")  | strpos(schnam, "DROP BACK IN")
 
 /*maybe*/
 drop if strpos(schnam, "BLIND") | strpos(schnam, "DEAF") | strpos(schnam, "SPECIAL NEEDS")
 drop if strpos(schnam, "CITRUS HEALTH") | strpos(schnam, "HIGHLAND PARK SIPP") | strpos(schnam, "GULF/LAKE ACADEMY") ///
 | strpos(schnam, "THE PORT ACADEMY") | strpos(schnam, "E SEAL") 
 drop if zip==32084  // this is the zip code with all the special ed schools





*******************************************************************************
****************some additional sample changes********************
***************************************************************************************

*drop virtual schools
drop if strpos(schnam, "VIRTUAL")  | strpos(schnam, "ONLINE") | strpos(schnam, "ADULT") 

replace stability = stability/100 if stability>1
replace FRL = FRL/100 if FRL>1

*claudia method variables;
g gradechange6= down60_any_mjrhwy_4*gradedum4 
g gradechange9=down60_any_mjrhwy_4*gradedum7 
g gradechange6or9=gradechange6
replace gradechange6or9=1 if gradechange9==1

tsset id year
sort id year


***********Identifying people who went downwind in 6th grade, went upwind in 6th grade, or stayed the same

local downwindvariable down60_any_mjrhwy_4 
*to downwind
	gen temp=0
	replace temp=1 if L.`downwindvariable'==0&`downwindvariable'==1&grade1==6&L.grade1==5
	egen todownwind=max(temp) if `downwindvariable'~=., by(id)  /*DES: without conditional on not missing then if you went to
	school outside of 0.4 miles but otherwise moved downwind: then you get a 1 to this value and are included in observations you are outside of 0.4 miles in the
	event study */
	label var todownwind "Moved from upwind to downwind school in 6th grade"
	drop temp
	gen group=.
	replace group=4 if todownwind==1

*to upwind
	gen temp=0
	replace temp=1 if L.`downwindvariable'==1&`downwindvariable'==0&grade1==6&L.grade1==5
	egen toupwind=max(temp) if `downwindvariable'~=., by(id) 
	label var toupwind "Moved from downwind to upwind school in 6th grade"
	drop temp
	replace group=3 if toupwind==1
	
*always upwind
	gen temp=0
	replace temp=1 if L.`downwindvariable'==0&`downwindvariable'==0&grade1==6&L.grade1==5
	egen allupwind=max(temp) if `downwindvariable'~=., by(id) 
	label var allupwind "In upwind school from 5th to 6th grade"
	drop temp
	replace group=1 if allupwind==1
	sum `downwindvariable' if allupwind==1
	
*always downwind
	gen temp=0
	replace temp=1 if L.`downwindvariable'==1&`downwindvariable'==1&grade1==6&L.grade1==5
	egen alldownwind=max(temp) if `downwindvariable'~=., by(id) 
	label var alldownwind "In doiwnwind school from 5th to 6th grade"
	drop temp
	replace group=2 if alldownwind==1
	sum `downwindvariable' if alldownwind==1
	
	label define definition 4 "To downwind" 3 "To upwind" 1 "Always upwind" 2 "Always downwind"
	label values group definition
	
	tab group `downwindvariable', missing
	
**identifying % years in school downwind before/after
	egen temp1=mean(`downwindvariable') if grade1==3|grade1==4|grade1==5, by(id)
	egen mean_downwindbefore=mean(temp1), by(id)
	egen temp2=mean(`downwindvariable') if grade1==6|grade1==7|grade1==8, by(id)
	egen mean_downwindafter=mean(temp2), by(id)  /*DES note: not exactly  # of years, if missing
	one or more grade, assigned the mean of the non missing values*/
	

drop temp*

***Identifying people who went downwind in 9th grade, went upwind in 9th grade, or stayed the same
*to downwind
	gen temp=0
	replace temp=1 if L.`downwindvariable'==0&`downwindvariable'==1&grade1==9&L.grade1==8
	egen htodownwind=max(temp) if `downwindvariable'~=., by(id) 
	label var htodownwind "Moved from upwind to downwind school in 9th grade"
	drop temp
	gen hgroup=.
	replace hgroup=4 if htodownwind==1
	
	
*to upwind
	gen temp=0
	replace temp=1 if L.`downwindvariable'==1&`downwindvariable'==0&grade1==9&L.grade1==8
	egen htoupwind=max(temp) if `downwindvariable'~=., by(id) 
	label var htoupwind "Moved from downwind to upwind school in 9th grade"
	drop temp
	replace hgroup=3 if htoupwind==1
	

	
*always upwind
	gen temp=0
	replace temp=1 if L.`downwindvariable'==0&`downwindvariable'==0&grade1==9&L.grade1==8
	egen hallupwind=max(temp) if `downwindvariable'~=., by(id) 
	label var hallupwind "In upwind school from 8th to 9th grade"
	drop temp
	replace hgroup=1 if hallupwind==1
	sum `downwindvariable' if hallupwind==1
	

	
*always downwind
	gen temp=0
	replace temp=1 if L.`downwindvariable'==1&`downwindvariable'==1&grade1==9&L.grade1==8
	egen halldownwind=max(temp) if `downwindvariable'~=., by(id) 
	label var halldownwind "In doiwnwind school from 8th to 9th grade"
	drop temp
	replace hgroup=2 if halldownwind==1
	sum `downwindvariable' if halldownwind==1
	

			
	
	label define hdefinition 4 "To downwind" 3 "To upwind" 1 "Always upwind" 2 "Always downwind"
	label values hgroup hdefinition
	

**identifying % years in school downwind before/after
	egen temp1=mean(`downwindvariable') if grade1==6|grade1==7|grade1==8, by(id)
	egen hmean_downwindbefore=mean(temp1), by(id)
	egen temp2=mean(`downwindvariable') if grade1==9|grade1==10|grade1==11, by(id)
	egen hmean_downwindafter=mean(temp2), by(id)  /*DES note: not exactly  # of years, if missing
	one or more grade, assigned the mean of the non missing values*/


	
	

	forv g=1/4 {
	forv i=1/10 {
		gen gradedum`i'_g`g'=(gradedum`i'==1&group==`g')
		*label var gradedum`i'_g`g' "Middle Mover groupXgrade"
		}
	}

	forv g=1/4 {
	forv i=6/10 {
		gen hgradedum`i'_g`g'=(gradedum`i'==1&hgroup==`g')
		label var hgradedum`i'_g`g' "High Mover groupXgrade"
		display "grade `i' group `g'"
		}
	}



*************************************************************
*Figure 2: Effects by % of Time Downwind
************************************************************


********now do for jenni method

preserve
drop if changein6or9id==0

	keep if (grade1 == 5 | grade1 == 6 | grade1 == 8 | grade1 == 9)

**By 10 ppt cutoffs:  
foreach var in  mjrhwy { //IS USH
forv i=10(10)100 {
	g downwindmore`i'`var'_4=0
	replace downwindmore`i'`var'_4=1 if dwind_any_win4_`var'>`i'/100 & dwind_any_win4_`var'!=.
	
	replace downwindmore`i'`var'_4=. if dwind_any_win4_`var'==. /*not necessary because schools more than half a mile have already been dropped I assume? */
			
	}
}

foreach var in mjrhwy { //  IS USH

replace downwindmore20`var'_4=0 if downwindmore30`var'_4==1 | downwindmore40`var'_4==1 | downwindmore50`var'_4==1 | downwindmore60`var'_4==1 | downwindmore70`var'_4==1
replace downwindmore30`var'_4=0 if downwindmore40`var'_4==1 | downwindmore50`var'_4==1 | downwindmore60`var'_4==1 | downwindmore70`var'_4==1
replace downwindmore40`var'_4=0 if downwindmore50`var'_4==1 | downwindmore60`var'_4==1 | downwindmore70`var'_4==1
replace downwindmore50`var'_4=0 if downwindmore60`var'_4==1 | downwindmore70`var'_4==1
replace downwindmore60`var'_4=0 if downwindmore70`var'_4==1
}

foreach var in mjrhwy{ //stdfcatread stdfcatmath IS USH 
estimates clear

gen filler=0

set more off

reghdfe avgfcat downwindmore20`var'_4 downwindmore30`var'_4 downwindmore40`var'_4 downwindmore50`var'_4 downwindmore60`var'_4 downwindmore70`var'_4 ///
 ${model3} gradedum4 gradedum7, absorb(zip grade1 year id) vce(cluster id zip)

estimates store downwind

coefplot (downwind), drop(_cons momedbyschool momblackbyschool mommarriedbyschool gradedum4 gradedum7 frl mover distmj* rdcnt* percenthispanic TeachDegree size stability) yline(0, lcolor(black))  pstyle(p1) lcolor(blue) ///
	rename(filler "<20%" downwindmore20`var'_4="20-30%" downwindmore30`var'_4="30-40%" downwindmore40`var'_4="40-50%" downwindmore50`var'_4="50-60%" downwindmore60`var'_4="60-70%" ///
	downwindmore70`var'_4=">70%" downwindmore80_5="80-100%" downwindmore90_5="90%") ytitle("Effect relative to <20% downwind") ///
	nokey vertical xtitle("Percent of Time Downwind") recast(connected) lcolor(gs0) mcolor(gs0)   ///
	ciopts(recast(rline) lpattern(dash_dot) lcolor(gs2) graphregion(color(white)) bgcolor(white))
}

graph save "$output/Figure2_schooltimedownwindjenni_`var'10percentbins70up$S_DATEy.gph", replace
graph export "$output/Figure2_schooltimedownwindjenni_`var'10percentbins70up$S_DATEy.tif", replace 

drop filler

restore
*/

******************************************************************************
*Figure 3: Event studies 6th grade
****************************************************************************

graph drop _all

*******************************
*** middle schools switchers limited, up/downwind movers:
*******************************
gen graphgrade=.
	label var graphgrade "Change in scores by grade"
	replace graphgrade=0 if grade1==5  

gen graphgrade_sel=.
	label var graphgrade_sel "Lower 95% CI"
	replace graphgrade_sel=0 if grade1==5

gen graphgrade_seu=.
	label var graphgrade_seu "Upper 95% CI"
	replace graphgrade_seu=0 if grade1==5

foreach out in avgfcat behaveincident rabsent { 

preserve
***Main model, no abnormal switchers (called "limited"):
keep if grade1>=2.5&grade1<=7.5
drop if group==.


drop  gradedum*_g1 gradedum*_g2
*can elimintate the following two lines of code to remove limited:
drop if mean_downwindbefore!=0&mean_downwindbefore!=1
drop if mean_downwindafter !=0&mean_downwindafter!=1


if "`out'"=="avgfcat" local title = "Average FCAT"
else if "`out'"=="stdfcatmath" local title = "Math Score"
else if "`out'"=="stdfcatread" local title = "Reading Score"
else if "`out'"=="behaveincident" local title = "Behavioral Incident 0/1"
else if "`out'"=="rabsent" local title = "Percent Days Absent"
else local title = "none"

estimates clear

display "**************************
display "`title'"
display"*****************************"
eststo: reghdfe `out'  gradedum1 gradedum2 gradedum3 gradedum4 gradedum5 gradedum6 gradedum1_g* gradedum2_g* gradedum4_g* gradedum5_g* gradedum6_g* $model3, absorb(zip year id) vce(cluster id zip)
esttab using "$output\clpeventstudyoutput`out'y_$S_DATE.rtf", ///
	compress replace b(4) se(4) onecell label nonumbers star(+ .1 * .05 ** .01 *** .001) ///
	title ("effect of `var' on outcomes") keep(gradedum1_g* gradedum2_g* gradedum4_g* gradedum5_g* gradedum6_g*)  
eststo clear 


forv g=3/4 {
		*generage coefficients to plot:
		foreach i in 1 2 4 5 6 {
			replace graphgrade=_b[gradedum`i'_g`g'] if grade1==`i'+2&group==`g'
			replace graphgrade_sel=_b[gradedum`i'_g`g']-1.96*_se[gradedum`i'_g`g'] if grade1==`i'+2&group==`g'
			replace graphgrade_seu=_b[gradedum`i'_g`g']+1.96*_se[gradedum`i'_g`g'] if grade1==`i'+2&group==`g'
		*move grade slightly y group so they don't overlap
			replace grade1=grade1-.1+`g'/20 if  grade1==`i'+2&group==`g'
		
		}
	}


collapse graphgrade graphgrade_sel graphgrade_seu, by(group grade1)


 twoway 	(connected graphgrade grade1 if group==3, 				 	sort lcolor(gs2) mcolor(gs2) lwidth(vthin) msymbol(D)) 	///
		(rcap graphgrade_sel graphgrade_seu grade1 if group==3, 	sort lcolor(gs2) lwidth(vthin) ) ///
		(connected graphgrade grade1 if group==4, 						sort lcolor(gs0) lwidth(vthin) mcolor(gs0)) ///
		(rcap graphgrade_sel graphgrade_seu grade1 if group==4, 	sort lcolor(gs0) lwidth(vthin) lpattern(dash_dot)), ///
		xtitle(Grade) xline(5.8, lcolor(black)) ///
		ytitle(Change from Grade 5) yline(0, lcolor(black)) ///
		title("`title'") ///  /*I can't figure out how to change the title depending on 1out' */
		graphregion(color(white)) ///
		legend(order( 1 "To upwind" 3 "To downwind") col(3)) ///
		name(mid`out'_lim, replace)
		restore
		}
	
	


		*graph save "$output\mid`out'_lim.gph", replace
		*graph export "$output\mid`out'_lim.pdf", replace
		

grc1leg  midavgfcat_lim midbehaveincident_lim midrabsent_lim, legendfrom(midavgfcat_lim) graphregion(color(white)) col(3)
		graph save "$output\Figure3_midcombined_limavgy.gph", replace
		graph export "$output\Figure3_midcombined_limavgy.pdf", replace		

	drop graphgrade graphgrade_sel graphgrade_seu
graph drop _all


*******************************
*** middle schools switchers limited, downwind movers only:
*******************************
gen graphgrade=.
	label var graphgrade "Change in scores by grade"
	replace graphgrade=0 if grade1==5  

gen graphgrade_sel=.
	label var graphgrade_sel "Lower 95% CI"
	replace graphgrade_sel=0 if grade1==5

gen graphgrade_seu=.
	label var graphgrade_seu "Upper 95% CI"
	replace graphgrade_seu=0 if grade1==5

foreach out in avgfcat behaveincident rabsent { 

preserve
***Main model, no abnormal switchers (called "limited"):
keep if grade1>=2.5&grade1<=7.5
drop if group==.|group==4 // drops downwind movers 


drop  gradedum*_g1 gradedum*_g2 gradedum*_g4 // rops downwind dummies
 *can elimintate the following two lines of code to remove limited:
drop if mean_downwindbefore!=0&mean_downwindbefore!=1
drop if mean_downwindafter !=0&mean_downwindafter!=1


if "`out'"=="avgfcat" local title = "Average FCAT"
else if "`out'"=="stdfcatmath" local title = "Math Score"
else if "`out'"=="stdfcatread" local title = "Reading Score"
else if "`out'"=="behaveincident" local title = "Behavioral Incident 0/1"
else if "`out'"=="rabsent" local title = "Percent Days Absent"
else local title = "none"

estimates clear

display "**************************
display "`title'"
display"*****************************"
eststo: reghdfe `out'  gradedum1 gradedum2 gradedum3 gradedum4 gradedum5 gradedum6 gradedum1_g* gradedum2_g* gradedum4_g* gradedum5_g* gradedum6_g* $model3, absorb(zip year id) vce(cluster id zip)
esttab using "$output\clpeventstudyoutput`out'_up_y_$S_DATE.rtf", ///
	compress replace b(4) se(4) onecell label nonumbers star(+ .1 * .05 ** .01 *** .001) ///
	title ("effect of `var' on outcomes") keep(gradedum1_g* gradedum2_g* gradedum4_g* gradedum5_g* gradedum6_g*)  
eststo clear 


forv g=3/3 {
		*generage coefficients to plot:
		foreach i in 1 2 4 5 6 {
			replace graphgrade=_b[gradedum`i'_g`g'] if grade1==`i'+2&group==`g'
			replace graphgrade_sel=_b[gradedum`i'_g`g']-1.96*_se[gradedum`i'_g`g'] if grade1==`i'+2&group==`g'
			replace graphgrade_seu=_b[gradedum`i'_g`g']+1.96*_se[gradedum`i'_g`g'] if grade1==`i'+2&group==`g'
		*move grade slightly y group so they don't overlap
			replace grade1=grade1-.1+`g'/20 if  grade1==`i'+2&group==`g'
		
		}
	}


collapse graphgrade graphgrade_sel graphgrade_seu, by(group grade1)


 twoway 	(connected graphgrade grade1 if group==3, 				 	sort lcolor(gs2) mcolor(gs2) lwidth(vthin) msymbol(D)) 	///
		(rcap graphgrade_sel graphgrade_seu grade1 if group==3, 	sort lcolor(gs2) lwidth(vthin) ) ///
		, ///
		xtitle(Grade) xline(5.8, lcolor(black)) ///
		ytitle(Change from Grade 5) yline(0, lcolor(black)) ///
		title("`title'") ///  /*I can't figure out how to change the title depending on 1out' */
		graphregion(color(white)) ///
		legend(order( 1 "To upwind") col(3)) ///
		name(mid`out'_lim, replace)
		restore
		}
	
		*(connected graphgrade grade1 if group==4, 						sort lcolor(navy) lwidth(vthin) mcolor(navy)) ///
		*(rcap graphgrade_sel graphgrade_seu grade1 if group==4, 	sort lcolor(navy) lwidth(vthin) lpattern(dash_dot)), ///
	


		*graph save "$output\mid`out'_lim.gph", replace
		*graph export "$output\mid`out'_lim.pdf", replace
		

grc1leg  midavgfcat_lim midbehaveincident_lim midrabsent_lim, legendfrom(midavgfcat_lim) graphregion(color(white)) col(3)
		graph save "$output\Figure3_midcombined_limavgy_up.gph", replace
		graph export "$output\Figure3_midcombined_limavgy_up.pdf", replace		

	drop graphgrade graphgrade_sel graphgrade_seu
graph drop _all

*******************************
*** middle schools switchers limited, downwind movers only:
*******************************
gen graphgrade=.
	label var graphgrade "Change in scores by grade"
	replace graphgrade=0 if grade1==5  /*DES: we need to drop one otherwise the group dummy is colinear with the student fixed effects,
	since students never change their group */

gen graphgrade_sel=.
	label var graphgrade_sel "Lower 95% CI"
	replace graphgrade_sel=0 if grade1==5

gen graphgrade_seu=.
	label var graphgrade_seu "Upper 95% CI"
	replace graphgrade_seu=0 if grade1==5

foreach out in avgfcat behaveincident rabsent { 

preserve
***Main model, no abnormal switchers (called "limited"):
keep if grade1>=2.5&grade1<=7.5
drop if group==.|group==3 // drops upwind movers 


drop  gradedum*_g1 gradedum*_g2 gradedum*_g3 // rops upwind dummies
 *can elimintate the following two lines of code to remove limited:
drop if mean_downwindbefore!=0&mean_downwindbefore!=1
drop if mean_downwindafter !=0&mean_downwindafter!=1


if "`out'"=="avgfcat" local title = "Average FCAT"
else if "`out'"=="stdfcatmath" local title = "Math Score"
else if "`out'"=="stdfcatread" local title = "Reading Score"
else if "`out'"=="behaveincident" local title = "Behavioral Incident 0/1"
else if "`out'"=="rabsent" local title = "Percent Days Absent"
else local title = "none"

estimates clear

display "**************************
display "`title'"
display"*****************************"
eststo: reghdfe `out'  gradedum1 gradedum2 gradedum3 gradedum4 gradedum5 gradedum6 gradedum1_g* gradedum2_g* gradedum4_g* gradedum5_g* gradedum6_g* $model3, absorb(zip year id) vce(cluster id zip)
esttab using "$output\clpeventstudyoutput`out'_down_y_$S_DATE.rtf", ///
	compress replace b(4) se(4) onecell label nonumbers star(+ .1 * .05 ** .01 *** .001) ///
	title ("effect of `var' on outcomes") keep(gradedum1_g* gradedum2_g* gradedum4_g* gradedum5_g* gradedum6_g*)  
eststo clear 


forv g=4/4 {
		*generage coefficients to plot:
		foreach i in 1 2 4 5 6 {
			replace graphgrade=_b[gradedum`i'_g`g'] if grade1==`i'+2&group==`g'
			replace graphgrade_sel=_b[gradedum`i'_g`g']-1.96*_se[gradedum`i'_g`g'] if grade1==`i'+2&group==`g'
			replace graphgrade_seu=_b[gradedum`i'_g`g']+1.96*_se[gradedum`i'_g`g'] if grade1==`i'+2&group==`g'
		*move grade slightly y group so they don't overlap
			replace grade1=grade1-.1+`g'/20 if  grade1==`i'+2&group==`g'
		
		}
	}


collapse graphgrade graphgrade_sel graphgrade_seu, by(group grade1)


 twoway (connected graphgrade grade1 if group==4, 						sort lcolor(gs0) lwidth(vthin) mcolor(gs0)) ///
		(rcap graphgrade_sel graphgrade_seu grade1 if group==4, 	sort lcolor(gs0) lwidth(vthin) lpattern(dash_dot)) ///
		, ///
		xtitle(Grade) xline(5.8, lcolor(black)) ///
		ytitle(Change from Grade 5) yline(0, lcolor(black)) ///
		title("`title'") ///  /*I can't figure out how to change the title depending on 1out' */
		graphregion(color(white)) ///
		legend(order( 1 "To downwind") col(3)) ///
		name(mid`out'_lim, replace)
		restore
		}
	

	


		*graph save "$output\mid`out'_lim.gph", replace
		*graph export "$output\mid`out'_lim.pdf", replace
		

*grc1leg midfcatNOCI_lim midavgfcat_lim midbehaveincident_lim midgraderep_lim, legendfrom(midavgfcat_lim) graphregion(color(white)) 
grc1leg  midavgfcat_lim midbehaveincident_lim midrabsent_lim, legendfrom(midavgfcat_lim) graphregion(color(white)) col(3)
		graph save "$output\Figure3_midcombined_limavgy_down.gph", replace
		graph export "$output\Figure3_midcombined_limavgy_down.pdf", replace		

	drop graphgrade graphgrade_sel graphgrade_seu
graph drop _all


***********************************************
****Doing it for a different reference year:
***********************************************

graph drop _all

*******************************
*** middle schools switchers limited:
*******************************
gen graphgrade=.
	label var graphgrade "Change in scores by grade"
	replace graphgrade=0 if grade1==3  /*DES: we need to drop one otherwise the group dummy is colinear with the student fixed effects,
	since students never change their group */

gen graphgrade_sel=.
	label var graphgrade_sel "Lower 95% CI"
	replace graphgrade_sel=0 if grade1==3

gen graphgrade_seu=.
	label var graphgrade_seu "Upper 95% CI"
	replace graphgrade_seu=0 if grade1==3

foreach out in avgfcat behaveincident rabsent { 

preserve
***Main model, no abnormal switchers (called "limited"):
keep if grade1>=2.5&grade1<=7.5
drop if group==.


drop  gradedum*_g1 gradedum*_g2
*can elimintate the following two lines of code to remove limited:
drop if mean_downwindbefore!=0&mean_downwindbefore!=1
drop if mean_downwindafter !=0&mean_downwindafter!=1


if "`out'"=="avgfcat" local title = "Average FCAT"
else if "`out'"=="stdfcatmath" local title = "Math Score"
else if "`out'"=="stdfcatread" local title = "Reading Score"
else if "`out'"=="behaveincident" local title = "Behavioral Incident 0/1"
else if "`out'"=="rabsent" local title = "Percent Days Absent"
else local title = "none"

estimates clear

display "**************************
display "`title'"
display"*****************************"
eststo: reghdfe `out'  gradedum1 gradedum2 gradedum3 gradedum4 gradedum5 gradedum6 gradedum2_g* gradedum3_g* gradedum4_g* gradedum5_g* gradedum6_g* $model3, absorb(zip year id) vce(cluster id zip)
esttab using "$output\clpeventstudyoutput`out'y_$S_DATE.rtf", ///
	compress replace b(4) se(4) onecell label nonumbers star(+ .1 * .05 ** .01 *** .001) ///
	title ("effect of `var' on outcomes") keep(gradedum2_g* gradedum3_g* gradedum4_g* gradedum5_g* gradedum6_g*)  
eststo clear 


forv g=3/4 {
		*generage coefficients to plot:
		foreach i in  2 3 4 5 6 {
			replace graphgrade=_b[gradedum`i'_g`g'] if grade1==`i'+2&group==`g'
			replace graphgrade_sel=_b[gradedum`i'_g`g']-1.96*_se[gradedum`i'_g`g'] if grade1==`i'+2&group==`g'
			replace graphgrade_seu=_b[gradedum`i'_g`g']+1.96*_se[gradedum`i'_g`g'] if grade1==`i'+2&group==`g'
		*move grade slightly y group so they don't overlap
			replace grade1=grade1-.1+`g'/20 if  grade1==`i'+2&group==`g'
		
		}
	}


collapse graphgrade graphgrade_sel graphgrade_seu, by(group grade1)


 twoway 	(connected graphgrade grade1 if group==3, 				 	sort lcolor(gs2) mcolor(gs2) lwidth(vthin) msymbol(D)) 	///
		(rcap graphgrade_sel graphgrade_seu grade1 if group==3, 	sort lcolor(gs2) lwidth(vthin) ) ///
		(connected graphgrade grade1 if group==4, 						sort lcolor(gs0) lwidth(vthin) mcolor(gs0)) ///
		(rcap graphgrade_sel graphgrade_seu grade1 if group==4, 	sort lcolor(gs0) lwidth(vthin) lpattern(dash_dot)), ///
		xtitle(Grade) xline(5.8, lcolor(black)) ///
		ytitle(Change from Grade 3) yline(0, lcolor(black)) ///
		title("`title'") ///  /*I can't figure out how to change the title depending on 1out' */
		graphregion(color(white)) ///
		legend(order( 1 "To upwind" 3 "To downwind") col(3)) ///
		name(mid`out'_lim, replace)
		restore
		}
	
	


		*graph save "$output\mid`out'_lim.gph", replace
		*graph export "$output\mid`out'_lim.pdf", replace
		

*grc1leg midfcatNOCI_lim midavgfcat_lim midbehaveincident_lim midgraderep_lim, legendfrom(midavgfcat_lim) graphregion(color(white)) 
grc1leg  midavgfcat_lim midbehaveincident_lim midrabsent_lim, legendfrom(midavgfcat_lim) graphregion(color(white)) col(3)
		graph save "$output\Figure3_midcombined_limavgy_alt.gph", replace
		graph export "$output\Figure3_midcombined_limavgy_alt.pdf", replace		

	drop graphgrade graphgrade_sel graphgrade_seu
graph drop _all





***************************************************************
*Figure 4: School Demographics over % downwind
********************************************************************
preserve

*tab FRL
*tab dwind_any_win4_mjrhwy

drop if changein6or9id==0

	keep if (grade1 == 5 | grade1 == 6 | grade1 == 8 | grade1 == 9)

*downwind bins;

gen dbins = 0.1 if (dwind_any_win4_mjrhwy >=0 & dwind_any_win4_mjrhwy <0.1) & dwind_any_win4_mjrhwy ~=.
replace dbins = 0.2 if (dwind_any_win4_mjrhwy >=0.1 & dwind_any_win4_mjrhwy <0.2) & dwind_any_win4_mjrhwy ~=.
replace dbins = 0.3 if (dwind_any_win4_mjrhwy >=0.2 & dwind_any_win4_mjrhwy <0.3) & dwind_any_win4_mjrhwy ~=.
replace dbins = 0.4 if (dwind_any_win4_mjrhwy >=0.3 & dwind_any_win4_mjrhwy <0.4) & dwind_any_win4_mjrhwy ~=.
replace dbins = 0.5 if (dwind_any_win4_mjrhwy >=0.4 & dwind_any_win4_mjrhwy <0.5) & dwind_any_win4_mjrhwy ~=.
replace dbins = 0.6 if (dwind_any_win4_mjrhwy >=0.5 & dwind_any_win4_mjrhwy <0.6) & dwind_any_win4_mjrhwy ~=.
replace dbins = 0.7 if (dwind_any_win4_mjrhwy >=0.6 & dwind_any_win4_mjrhwy <0.7) & dwind_any_win4_mjrhwy ~=.
replace dbins = 0.8 if (dwind_any_win4_mjrhwy >=0.7 & dwind_any_win4_mjrhwy <0.8) & dwind_any_win4_mjrhwy ~=.
replace dbins = 0.9 if (dwind_any_win4_mjrhwy >=0.8 & dwind_any_win4_mjrhwy <0.9) & dwind_any_win4_mjrhwy ~=.
replace dbins = 1.0 if (dwind_any_win4_mjrhwy >=0.9 & dwind_any_win4_mjrhwy <=1.0) & dwind_any_win4_mjrhwy ~=.

tab dbins, missing

sum dbins dwind_any_win4_mjrhwy


collapse FRL momblackbyschool mommarried momedbyschool (sum) weight=id, by(dbins) fast

rename dbins dwind_any_win4_mjrhwy


*linear fit

# delimit ;
	twoway (scatter FRL dwind_any_win4_mjrhwy, mcolor(gs0)) (lfitci FRL dwind_any_win4_mjrhwy [aw=weight], ciplot(rcap) lcolor(gs10)), 		 
			 ytitle("Mean Characteristic", size(small)) xtitle(" ")
			 title("Percent FRL Across Percent Downwind" , size(medsmall) color(gs0))
			 yscale(range(0))  ylabel(0(.2)1)
			 saving("$output/percent frl across percent downwind mjrhwy.gph", replace)
			 graphregion(color(white))
			 yscale(range(0)) legend(off)
			 ;
	# delimit cr

# delimit ;
	twoway  (scatter momedbyschool dwind_any_win4_mjrhwy, mcolor(gs0)) (lfitci momedbyschool dwind_any_win4_mjrhwy, ciplot(rcap) lcolor(gs10) ),			 
			 yscale(range(0))   ylabel(0(3)18) ytitle(" ") xtitle(" ")
			 title("Years of Maternal Education Across Percent Downwind" , size(medsmall) color(gs0))
			 saving("$output/mom ed across percent downwind mjrhwy.gph", replace)
			 graphregion(color(white))
			 yscale(range(0))  legend(off)
			 ;
	# delimit cr


# delimit ;
	twoway (scatter momblackbyschool dwind_any_win4_mjrhwy, mcolor(gs0)) (lfitci momblackbyschool dwind_any_win4_mjrhwy, ciplot(rcap) lcolor(gs10) ) , 			 
			
			 ytitle("Mean Characteristic", size(small)) xtitle("Percent Downwind")
			 yscale(range(0))   ylabel(0(.2)1) 
			 title("Percent Black Across Percent Downwind" , size(medsmall) color(gs0))
			 saving("$output/percent black across percent downwind mjrhwy.gph", replace)
			 graphregion(color(white))
			 yscale(range(0))  legend(off)
			 ;
	# delimit cr

# delimit ;
	twoway (scatter mommarried dwind_any_win4_mjrhwy, mcolor(gs0)) (lfitci mommarried dwind_any_win4_mjrhwy, ciplot(rcap) lcolor(gs10) )	,	 
			  yscale(range(0))   ylabel(0(.2)1)  ytitle(" ") xtitle("Percent Downwind")
			 title("Percent in Married families Across Percent Downwind" , size(medsmall) color(gs0))
			 saving("$output/percent married moms percent downwind mjrhwy.gph", replace)
			 graphregion(color(white)) legend(off)
			 yscale(range(0))  
			 ;
	# delimit cr

# delimit ;
	graph combine "$output/percent frl across percent downwind mjrhwy.gph" "$output/mom ed across percent downwind mjrhwy.gph"
	"$output/percent black across percent downwind mjrhwy.gph"  "$output/percent married moms percent downwind mjrhwy.gph" ,
		saving("$output/Figure4.gph", replace) 
		graphregion(color(white))
		;

	# delimit cr



restore



*************************************************************
*Figure 5: Effects by Traffic
*************************************************************
preserve



drop if changein6or9id==0

	keep if (grade1 == 5 | grade1 == 6 | grade1 == 8 | grade1 == 9)


local downwindvariable closetohwydwind60_4  //down60_any_mjrhwy_4 down60_any_aadt25k_4
summarize AADT1_mjrhwy if `downwindvariable'!=., meanonly

replace closetohwydwind60_4=. if mitohwy>0.4 //eliminating the >0.4 mi comparison group makes traffic density more important
//down60_any_aadt48k??

gen aadtless25=0
replace aadtless25 = 1 if AADT1_mjrhwy<25000 

gen aadt25to50 =0
replace aadt25to50 = 1 if AADT1_mjrhwy>=25000 & AADT1_mjrhwy<50000 

gen aadt50to75 =0
replace aadt50to75 = 1 if AADT1_mjrhwy>=50000 & AADT1_mjrhwy<75000 

gen aadt75more = 0
replace aadt75more = 1 if AADT1_mjrhwy>=75000


gen aadtless50=0
replace aadtless50 = 1 if AADT1_mjrhwy<50000


foreach i in aadtless25 aadt25to50 aadt50to75 aadt75more aadtless50 { // addt75to100 addt100more addt125to150 addt150to175 addt175to200 addt200more 
g `i'i=`i'*`downwindvariable'
}


tab aadtless50i
tab aadt50to75i
tab aadt75morei

set more off

estimates clear
reghdfe avgfcat aadtless50i aadt50to75i aadt75morei aadt50to75 aadt75more ${model3}, absorb(zip grade1 year id) vce(cluster id zip) //aadt25to50i aadt75to100i addt100morei aadt25to50i addt100to125i addt125to150i addt150to175i addt175to200i addt200morei
estimates store aadt

*local reference = _b[aadtless25]

coefplot (aadt),  drop(_cons momedbyschool momblackbyschool mommarriedbyschool frl mover aadt25to50 aadt50to75 aadt75more distmj* rdcnt* percenthispanic TeachDegree size stability) yline(0, lcolor(gs10)) ///
 pstyle(p1) lcolor(gs0) mcolor(gs0) ///
rename(aadtless50i="AADT<50,000"  aadt50to75i="AADT 50,000-75,000" aadt75morei="AADT>75,000") ytitle("Average FCAT") ///
nokey vertical xtitle("Traffic") recast(connected) lcolor(gs0)  lwidth(vthin) coeflabels(aadtless50="AADT<50,000" aadt25to50i="AADT 25,000-50,000" aadt50to75i="AADT 50,000-75,000" aadt75morei="AADT>75,000")  ///
ciopts(recast(rline)  lpattern(dash_dot) lcolor(gs2) lwidth(vthin) graphregion(color(white)) bgcolor(white))
		
		graph save "$output\Figure5_traffic_coefploty_$S_DATE.gph", replace
		graph export "$output\Figure5_traffic_coefploty_$S_DATE.pdf", replace
restore

*/

*************************************************************
* Appendix Figure A1: Event studies All movers
************************************************************
preserve

********************************************************************************
scalar downspell = 2	// Indicates the min number of periods a student must stay in downwind schools once she switches
scalar upspell = 3	// Indicates the min number of periods a student must have spent in upwind schools before switching to downwind 
*including the -1 reference period
scalar lag = 2			// This indicates the number of leads and lags in the event study (not counting reference period)

global model3 momedbyschool momblackbyschool mommarriedbyschool percenthispanic frl mover TeachDegree size stability distmj* rdcnt_win10_mjrhwy 
********************************************************************************




order id year school 
qui sum id
scalar maxid = r(max)


drop if down60_any_mjrhwy_4==.

cap drop yearnew
bysort id (year): gen yearnew = _n
tab year yearnew, missing

sort id yearnew
xtset id yearnew
tsspell down60_any_mjrhwy_4 

bysort id _spell: egen length = max(_seq)	// gives the length of each spell


bysort id (year): gen switch = down60_any_mjrhwy_4 - down60_any_mjrhwy_4[_n-1]
bysort id (year): replace switch= 0 if _n==1  // don't count the first year as a switch


replace switch=0 if switch==.



gen yearswitch = yearnew if down60_any_mjrhwy_4==1 & length>=downspell & switch==1 & length[_n-1]>=upspell // last condition indicates upwind spell length


egen ctswitch = count(yearswitch), by(id)
tab ctswitch, missing


 egen trtswitch = min(switch) if ctswitch==0, by(id)
 
 tab trtswitch, missing
 
  
 drop if trtswitch==-1  


  egen parswitch = max(switch) if ctswitch==0, by(id)
 *this drops those who switch from upwind to downwind, but don't have a qualifying event

 tab parswitch, missing
 
 drop if parswitch==1

assert switch==0 if ctswitch==0


*note if switch==1 mover equals 1 but not always the other way

assert mover==1 if switch==1


	
*drop those with multiple qualifying events for now;  

drop if ctswitch>1


egen eventyear = max(yearswitch), by(id)

gen timeto = yearnew-eventyear

assert timeto==. if eventyear==.

order id yearnew eventyear timeto down60_any_mjrhwy_4 switch yearswitch ctswitch 



gen pred1 = timeto==-1 
scalar s = `=lag'+1
forvalues j=1/`=s' {
	local k=`j'+1
	local f=`j'-1
	gen pred`k'= timeto==-1*`k' 
	gen postd`f'=timeto==`f'
					}
sum pred* postd*

					
local z = `=lag'+2
replace pred`z'=1 if timeto<-1*(`z') & eventyear~=.
replace postd`=lag'=1 if timeto>`=lag' & eventyear~=.

rename pred`z' farpre
rename postd`=lag' farpost

*making endpoints inclusive
replace pred3=1 if farpre ==1
replace postd1=1 if farpost==1

sum farpre pred* postd* farpost
* these all check out

************************************
// EVENT STUDY
************************************
// regress Y on leads and lags, with student by year, school and zipcode fixed effects
//egen studyear = group(id year)
set more off
estimates clear
set emptycells drop

// The following code runs the event study, creates a dataset with the estimates, 
// sets the estimate of the reference year to 0, and cretes id variables used in the graph
// If changing the number of lags, they have to be coded manually right now

*make year -1 be the reference year;

*drop pred1

gen ref = 1


*do average fcat no CI first.  




foreach vv in avgfcat behaveincident rabsent {

if "`vv'"=="avgfcat" local title = "Average FCAT"
else if "`vv'"=="behaveincident" local title = "Behavioral Incident 0/1"
else if "`vv'"=="rabsent" local title = "Percent Days Absent"
else local title = "none"

display " downwind regression on `vv', Event sample, "
estimates clear
	
display "Eventstudy on `vv'"

reghdfe `vv' pred3 pred2 ref postd0 postd1 $model3, absorb(zip grade1 year id) vce(cluster id zip) 


coefplot, recast(connected) name(`vv', replace) lcolor(gs0) mcolor(gs0) lwidth(vthin) ///
	legend(order(2 "Effect of Moving Downwind")) title("`title'", size(medsmall)) ///
	ciopts(recast(rcap) lpattern(dash_dot) lcolor(gs2) lwidth(vthin)) ///
	graphregion(color(white)) bgcolor(white) ///
	vertical coeflabels(pred3="-3" pred2="-2" ref="-1" postd0="0" postd1="+1") keep(pred* ref postd*) ///
	omitted xline(3.6, lcolor(gs10)) yline(0, lwidth(medthin) lcolor(gs10))  xtitle("Event Year")
}
grc1leg avgfcat behaveincident rabsent, graphregion(color(white)) col(3)
		graph save "$output\FigureA2_combineallmv8-9-18y.gph", replace
		graph export "$output\FigureA2_combineallmv8-9-18y.pdf", replace	


restore


/*

*************************************************************
*Bonus table: robustness of core results to different clustering 
*************************************************************


eststo clear 

* Jenni Method


drop if changein6or9id==0

	keep if (grade1 == 5 | grade1 == 6 | grade1 == 8 | grade1 == 9)

foreach var in down60_any_mjrhwy_4  dwintensity4h  { 
	foreach i in  avgfcat behaveincident rabsent {  

	sum `i' if down60_any_mjrhwy_4~=. & grade1~=. & size~=. & zip~=.
	local m = r(mean)
	
eststo: reghdfe `i' `var' ${model3} gradedum4 gradedum7, absorb(zip year id) vce(cluster school)  
	estadd scalar mean_outcome=`m'


eststo: reghdfe `i' `var' ${model3} gradedum4 gradedum7, absorb(zip year id) vce(cluster id)  
	estadd scalar mean_outcome=`m'

eststo: reghdfe `i' `var' ${model3} gradedum4 gradedum7, absorb(zip year id) vce(cluster school id)  
	estadd scalar mean_outcome=`m'


eststo: reghdfe `i' `var' ${model3} gradedum4 gradedum7, absorb(zip year id) vce(cluster zip )  
	estadd scalar mean_outcome=`m'


eststo: reghdfe `i' `var' ${model3} gradedum4 gradedum7, absorb(zip year id) vce(cluster zip id)  
	estadd scalar mean_outcome=`m'


}


  esttab using "$output/cluster_`var'_$S_DATE.rtf", ///
	compress replace b(4) se(4) onecell label nonumbers star(* .1 ** .05 *** .01) ///
	title ("effect of `var' on outcomes") keep(`var') scalars(mean_outcome)
	
eststo clear 
}




log close

stop
